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1.0  ABSTRACT 


This  report  sunmarizes  results  obtained  tinder  contract  Bonr  5054(00). 

The  objective  of  this  work  was  the  construction  of  a  computer  program  for  the 
solution  of  the  problem  of  an  arbitrary  three-dimensional  body  with  zero 
translational  velocity  performing  small  steady  harmonic  oscillations  in  the 
presence  of  an  otherwise  undisturbed  free  surface,  which  bounds  the  fluid  in 
which  the  body  is  immersed.  Because  the  sponsorship  was  terminated,  this  work 
was  not  carried  to  completion.  The  results  that  were  obtained  are  contained 
in  this  report. 

The  solution  of  the  above  problem  consists  of  two  principal  tasks:  the 
development  of  certain  analytic  formulas  and  the  incorporation  of  these  formulas 
into  a  working  computer  program.  Almost  all  of  the  necessary  formulas  were 
derived  and  tested  numerically,  but  they  could  not  be  incorporated  into  a  pro¬ 
gram  in  the  time  available.  The  subsequent  sections  of  this  report  describe 
the  general  scheme  for  solving  the  problem  of  interest,  present  the  formulas 
that  were  derived  to  implement  this  scheme,  and  point  out  areas  where  additional 
formulas  need  to  be  developed. 
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5-0  PRINCIPAL  NOTATION 


A 

A 

o 


a,b,c 


mn 


■kW 

B 


b 

mn 

bk(P) 

Ck(B) 

Cfc(p) 

C«m 

mn 


Ei 

f 

fk(0 

o(a,p) 


equals  vR,  the  dimensionless  horizontal  distance  between  a  point 
source  and  a  field  point  where  the  potential  is  to  be  evaluated 

equals  vRq,  the  dimensionless  horizontal  distance  between  a  field 
point  and  the  centroid  of  an  element 

normal  velocity  at  the  control  point  of  the  i-th  surface  element  due 
to  a  unit  source  density  on  the  j-th  element 

coordinates  of  the  location  of  a  point  source 

m  =  1,2,3,  n  =  1,2,3;  components  of  the  unit  vectors  along  the  axes 
of  the  coordinate  system  of  an  element,  equation  (70) 

coefficient  functions  defined  by  equation  (58) 

equals  vh,  the  dimensionless  vertical  distance  between  the  image 
of  a  point  source  and  a  field  point  where  the  potential  is  to  be 
evaluated 

equals  vhQ,  the  dimensionless  vertical  distance  between  a  field 
point  and  the  image  of  the  centroid  of  an  element 

with  superscripts  x,  y,  and  z,  these  are  coefficients  in  the 
multipole  expansion 

coefficient  functions  defined  by  equation  (63) 

coefficient  functions  defined  by  equations  (48)  and  (52) 

coefficient  functions  defined  by  equation  (66) 

with  superscripts  x,  y,  and  z,  these  are  coefficients  in  the 
multipole  expansion 

exponential  integral,  equation  (40) 

function  that  specifies  normal  velocity  distribution  on  the  body 
surface 

coefficient  functions  defined  by  equations  (64)  and  (67) 

leading  term  in  the  expansion  of  section  7*5  for  the  potential, 
equation  (68) 
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g(A,p)  function  defined  in  terms  of  G(A,p)  by  equation  (69) 

H  weights  in  Laguerre-Gauss  quadrature 

Hankel  function  of  first  kind  of  order  zero 
o 

H(A, p)  leading  term  in  the  expansion  of  section  7*5  for  the  A-derivatlve 
of  the  potential, equation  (68) 

h  vertical  distance  between  the  image  of  a  point  source  in  the  free 

surface  and  a  field  point  where  the  potential  is  to  be  evaluated, 
figure  2 

hQ  vertical  distance  between  a  field  point  and  the  image  of  the  centroid 

of  an  element 

h(A,p)  function  defined  in  terms  of  H(A, p)  by  equation  (69) 

I  ,J  Bessel  functions  of  first  and  second  kinds  of  order  zero 

o’  o 

I  normalized  moments  of  the  area  of  a  quadrilateral  element  about  its 

mn 

centroid,  equation  (72) 

i, j  subscripts  denoting  quantities  associated  with  the  i-th  and  j-th 

elements,  respectively  (other  uses  of  these  subscripts  with  certain 
variables  are  detailed  explicitly  in  this  section) 

i,j,k  unit  vectors  along  the  axes  of  the  x,y,z-coordinate  system 

k  integer  subscript  used  in  several  ways,  especially  to  denote  terms 

of  expansions 

M  integer  denoting  the  order  of  an  expansion,  usually  the  upper  limit 

of  a  summation 

N  number  of  elements  used  to  approximate  the  body  surface 

n  distance  normal  to  the  body  surface 

tm h 

n  unit  normal  vector  to  the  body  surface 

P  point  off  the  body  surface  where  the  potential  is  to  be  evaluated 

p  point  on  the  body  surface  where  the  potential  is  to  be  evaluated; 

with  subscript  i,  the  control  point  of  the  i-th  element 

q  point  where  a  source  is  located,  especially  a  point  on  the  body 

surface 
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horizontal  distance  between  a  point  source  and  a  field  point  where 
the  potential  is  to  be  evaluated,  figure  2 

horizontal  distance  between  a  field  point  and  the  centroid  of  an 
element 

distance  between  a  point  source  and  a  field  point  where  the  potential 
is  to  be  evaluated,  figure  2 

distance  between  the  image  of  a  point  source  in  the  free  surface  and 
a  field  point  where  the  potential  is  to  be  evaluated,  figure  2 

mean  position  of  the  surface  of  the  body  that  is  performing  oscillations 

principal  use  is  as  the  maximum  dimension  of  a  quadrilateral  element) 
in  certain  sections  also  used  as  time,  as  integration  variable,  and 
with  subscript  j  as  abscissas  for  Laguerre -Gauss  quadrature 

with  subscripts  x,y,z,  denotes  velocity  component 

velocity  at  the  control  point  of  the  i-th  surface  element  due  to  a 
unit  source  density  on  the  j-th  element 

vector  velocity 

Cartesian  coordinates.  The  mean  position  of  the  free  surface  is  the 
plane  y  =  0,  and  'the  flow  field  is  the  half-space  y  <  0.  Used 
to  denote  coordinates  of  a  field  point  where  the  potential  is  to  be 
evaluated 

coordinates  of  the  centroid  of  a  quadrilateral  element 
area  of  a  quadrilateral  element 
B  -  pA 

<^/g>  where  g  is  acceleration  of  gravity.  This  is  the  spatial 
circular  frequency  of  the  motion 

coordinate  system  based  on  a  quadrilateral  element 

a  particular  value  of  B/A  about  which  the  potential  function  is 
expanded 

fluid  motion  is  a  pure  harmonic  in  time  with  circular  frequency 
equal  to  a.  Also  used  as  surface  source  density. 


<b 


<p 


% 

<pl 

% 

% 

V2 


integrated  potential  of  a  quadrilateral  element.  Also  time- 
dependent  potential 

potential  at  the  control  point  of  the  i-th  element  due  to  a  unit 
source  density  on  the  j-th  element 

velocity  potential  after  harmonic  time  dependence  has  been  removed. 
Used  in  various  places  to  denote  various  potentials 

potential  defined  by  equation  (27) 

potential  defined  by  equation  (30) 

potential  defined  by  equation  (26) 

potential  defined  by  equation  (29) 

potential  of  an  oscillating  point  source 

Laplacian  operator 
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6.0  THE  METHOD  OF  SOLUTION 


6.1  General  Description  of  the  Method 

For  the  past  several  years  there  has  been  a  continuous  effort  at  this 
facility  directed  towards  the  solution  of  certain  three-dimensional  problems 
of  fluid  dynamics.  The  first  problem  studied  was  that  of  the  potential  flow 
about  an  arbitrary  three-dimensional  body  in  an  infinite  fluid.  This  work 
was  successfully  completed,  and  a  description  of  the  method  of  solution  is 
contained  in  references  [l]  and  [2].  A  subsequent  investigation  concerned  the 
solution  of  the  Helmholtz  equation  governing  the  scattering  and  radiation  of 
acoustic  waves  by  an  aribtrary  three-dimensional  body.  A  description  of  the 
method  of  solution  is  contained  in  reference  [3]. 

The  two  above-mentioned  problems  are  logically  similar,  and  the  same 
general  approach  was  adopted  to  effect  each  of  their  solutions.  This  same 
approach  can  be  used  to  solve  the  problem  of  interest  here.  In  each  case  it 
is  desired  to  solve  for  a  scalar  velocity  potential  that  satisfies  a  partic¬ 
ular  partial  differential  equation  together  with  certain  auxiliary  conditions 
and  whose  normal  derivative  is  specified  on  the  body  surface.  The  solution 
method  is  based  on  a  point  source  function  that  identically  satisfies  the 
partial  differential  equation  and  the  auxiliary  conditions.  A  distribution 
of  the  appropriate  type  of  source  density  is  assumed  to  lie  on  the  body  sur¬ 
face.  Applying  the  normal  derivative  condition  on  the  body  surface  then  leads 
to  a  Fredholm  integral  equation  of  the  second  kind  for  the  unknown  source 
density  distribution.  The  only  difference  between  the  three  above-mentioned 
problems  lies  in  the  nature  of  the  three  different  point-source  functions  that 
are  appropriate  to  the  various  problems. 

The  approximate  solution  of  the  integral  equation  is  effected  in  the 
following  manner.  The  body  surface  is  approximated  by  a  large  number  of  small 
plane  surface  elements,  over  each  of  which  the  source  density  is  assumed 
constant.  On  each  element  a  control  point  is  selected  at  which  the  boundary 
condition  is  to  be  satisfied.  The  basic  formulas  for  the  point  source  potential 
and  its  gradient  are  integrated  over  each  element  to  obtain  the  effects  of  the 
elements  at  each  others '  control  points  per  unit  value  of  source  density.  In 
particular,  a  matrix  is  obtained  whose  entries  are  the  normal  derivatives  of 
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potential  induced  by  the  elements  at  each  others'  control  points  for  a  unit 
value  of  source  density.  This  is  the  coefficient  matrix  for  a  set  of  linear 
algebraic  equations  for  the  unknown  values  of  the  source  density  on  the  ele¬ 
ments.  The  right-hand  sides  of  these  equations  are  the  prescribed  values  of 
the  normal  derivative  of  the  potential  at  the  control  points.  Once  this  set 
of  equations  has  been  solved  for  the  values  of  the  source  density  on  the 
elements,  all  other  quantities  of  interest  can  be  calculated  fairly  easily  at 
any  desired  point.  The  key  operation  in  this  solution  scheme  is  the  integra¬ 
tion  of  the  basic  formulas  for  the  basic  point  source  potential  and  its  gradi¬ 
ent  over  a  plane  surface  element.  This  is  also  the  only  operation  that  is 
different  for  the  three  above-mentioned  problems.  Thus  the  solution  of  the 
problem  of  oscillatory  motion  near  a  free  surface  is  primarily  a  matter  of 
performing  this  integration  numerically. 

6.2  Mathematical  Statement  of  the  Problem 

The  problem  of  a  body  with  zero  translational  velocity  performing  small 
steady  oscillations  in  the  presence  of  a  free  surface  is  well  known.  The 
mathematical  formulation  is  given  in  references  [4]  and  [5L  and  results  will 
simply  be  stated  here  with  no  attempt  at  a  derivation. 

Cartesian  coordinates  x,y,  z  are  assumed,  and  the  undisturbed  location 
of  the  free  surface  is  taken  as  the  plane  y  =  0.  The  fluid  motion  takes 
place  in  the  half  space  y  <  0,  (figure  l).  In  discussing  this  problem 
distance  perpendicular  to  the  free  surface,  i.e.,  y  distance,  is  denoted 
vertical  distance  and  distance  parallel  to  the  free  surface,  i.e.,  x  and 
z  distance,  is  denoted  horizontal  distance.  The  surface  of  the  body  that  is 
performing  small  oscillations  is  denoted  S.  Specifically,  S  represents 
the  mean  position  of  the  body  surface,  and  the  oscillations  take  place  about 
this  surface.  In  particular,  the  surface  S  is  independent  of  time.  The 
body  represented  by  S  may  be  surface-piercing,  as  shown  in  figure  1,  or 
completely  submerged.  Moreover,  S  may  be  multiply-connected  and  represent 
an  ensemble  of  bodies. 

The  fluid  velocity  field  in  the  half-space  y  <  0  is  assumed  to  be 
irrotational.  Thus,  the  velocity  is  equal  to  the  negative  gradient  of  a  scalar 
potential  function  <t(x,y,z,t),  which  is  a  function  of  the  time  t  as  well 
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Figure  1.  -  A  three-dimensional  body  near  a  tree-surface. 


as  of  position.  By  assumption  the  fluid  motion  is  harmonic  with  a  single 
circular  frequency  a.  Thus,  the  time -dependent  potential  <j(x,y,z,t)  can 
be  written 

4(x,y,z,t)  =  Re[cp(x,y,z)e"10t]  (l) 

where  the  potential  cp(x,y,z)  is  independent  of  time  and  is  thus  denoted  the 
steady  potential.  It  is  the  potential  cp  that  must  be  calculated. 


According  to  references  [4]  and  [5]  cp  must  satisfy  the  following 
equations : 


V2  cp  =  0 

for 

y  <  0 

(2) 

0 

11 

0- 

> 

1 

for 

y  =  0 

(3) 

(grad  cp)  =  0 

00 
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2  2  .  2 
p  =  x  +  z 


^  =  f  (S)  on  S 


The  parameter  v  is  defined  as 


where  g  is  the  acceleration  of  gravity.  In  (6)  a  denotes  distance  normal 


to  S. 


iHe  physical  significance  of  the  above  equations  is  as  follows.  Equation 

(2)  is  the  partial  differential  equation  for  9  and  expresses  the  conditions 
of  incompressibility  and  irrotationality.  Equations  (3),  (4),  and  (5)  are  the 
auxiliary  conditions  on  9.  Equation  (3)  is  the  linearized  free-surface 
condition,  which  requires  that  the  pressure  on  the  free-surface  be  constant. 
Equation  (4)  expresses  the  vanishing  of  the  disturbance  at  infinite  depth, 
and  equation  (5)  is  the  radiation  condition  that  requires  the  disturbance 

to  be  an  outgoing  wave  at  infinite  horizontal  distance.  Equation  (6)  is  the 
boundary  condition  on  S.  It  expresses  the  fact  that  the  normal  fluid  velocity 
on  S  must  be  specified  as  a  function  f(S)  of  position  on  the  surface. 

Often  the  fluid  normal  velocity  is  specified  as  equal  to  the  normal  velocity 
of  the  surface  S.  However,  in  the  case  of  a  known  incident  wave,  9  denotes 
the  disturbance  potential  due  to  the  body,  and  the  boundary  condition  expresses 
the  fact  that  the  normal  velocity  of  the  disturbance  must  cancel  that  of  the 
incident  wave  on  the  body  surface. 

6.3  The  Oscillating  Point-Source  Potential 

The  method  of  solving  the  mathematical  problem  defined  by  equations  (2), 

(3) ,  (*0,  (5)>  and  (6)  is  based  on  a  so-called  elementary  solution  or  point- 
source  solution.  The  point-source  potential  is  defined  as  the  one  that 
identically  satisfies  the  auxiliary  conditions  (3),  (4),  and  (5),  and  that 
satisfies  the  homogeneous  partial  differential  equation  (2)  throughout  space 
except  at  one  point,  where  it  is  singular.  The  singularity  is  such  that  the 
point-source  potential  satisfies  equation  (2)  with  the  right  side  replaced  by 


the  Dirac  delta  function  (a  function  zero  everywhere  except  at  one  point  where 
it  is  infinite  and  such  that  the  integral  of  the  function  over  any  volume  that 
includes  the  singularity  is  unity).  The  point  where  the  point-source  potential 
is  singular  is  the  location  of  the  point  source. 


Expressions  for  the  point-source  potential  are  given  in  reference  [4], 
If  the  source  is  located  at  the  point  (a,b,c)  its  potential  at  a  general 
point  (x,y,z)  is 


<p  =  1  +  1_+  2vevy  f  dy  +  2irivev^yfb^  H^(vR)  (8) 

Ts  r  r^  J  ri  o  '  ' 

00 

where 

r2  =  (x  -  a)2  +  (y  -  b)2  +  (z  -  c)2  (9) 

r2  =  (x  -  a)2  +  (y  +  b)2  +  (z  -  c)2  (10) 

K2  *  (x  -  a)2  +  (z  -  c)2  (11) 

and  where  =  J  +  IY  is  the  Hankel  function  of  the  first  kind  of  order 

o  00 

zero.  (J  and  Y  are  the  ordinary  Bessel  functions  of  order  zero. )  The 
00 

physical  significances  of  the  three  distances  r,  r^,  and  R  are  evident 
from  equations  (9),  (10),  and  (ll),  and  they  are  illustrated  in  figure  2.  The 
distance  r  is  the  distance  between  the  point  source  and  the  field  point  where 
the  potential  is  evaluated.  Similarly,  r^  is  the  distance  from  the  field 
point  to  the  image  point  of  the  source  in  the  free  surface.  Finally,  R  is 
the  horizontal  distance  between  source  and  field  point.  The  point  source 
whose  potential  is  given  by  (8)  is  defined  as  a  source  of  unit  strength. 


For  ordinary  problems  governed  by  Laplace's  equation  with  no  free  surface 
the  point-source  potential  is  simply  l/r.  Hie  first  two  terms  of  (8)  are  seen 
to  be  the  potentials  of  l/r- type  point  sources  —  one  located  at  the  point 
source  and  one  at  its  image  in  the  free  surface.  The  integral  term  in  (8)  may 
be  written 


2veVy 


/» 

■L  A 


-vu 


du 


\]  R2  +  (u+y+b)2 


(12) 
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Figure  2.  •  Distances  between  a  field  point  where  potential  is  evaluated  and 
a  point  source  and  its  image. 

The  square  root  on  the  right  side  of  (12)  is  the  distance  between  the  point 
(x,y,z)  and  the  point  (a,-b-u,  c).  Thus,  this  integral  represents  the 
potential  of  a  l/r-type  line  source  of  exponentially  decaying  strength  start¬ 
ing  at  the  image  point  (a,  -b,  c)  of  the  point  source  and  running  vertically 
downward  throu^i  the  free  surface  to  y  = 

For  large  values  of  vR,  the  function  H^(vR)  oscillates  with 
increasing  horizontal  distance  R  at  a  circular  frequency  of  v.  Bius  v 
is  denoted  the  spatial  circular  frequency  of  the  motion,  and  its  relation  to 
the  temporal  frequency  a  is  given  by  (7)* 


6.4  Representation  of  the  Solution  by  a  Source  Distribution 
on  the  Body  Surface 


The  unit  point-source  potential,  which  for  the.  present  problem  is  given 
by  (8),  is  a.  function  of  the  field  point  P  with  coordinates  x,y,z  where 
the  potential  is  evaluated  and  of  the  point  q  with  coordinates  a,b,c  where 
th  source  is  located.  Accordingly,  the  potential  of  a  unit  point  source  may 
be  written  cpg(P,q).  Now  let  the  source  point  q  be  restricted  to  lie  on  the 
surface  S.  (Figure  3.)  A  distribution  of  source  density  of  strength  cx(q) 

P 

f 
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Figure  3.  -  Illustration  of  quantities  appearing  in  the  integral  equation. 


is  assumed  to  lie  on  the  surface  ?.  The  potential  at  P  due  to  this 
distribution  is 


ff(p)  =^cPs(I,>q)°((l)dS  (13) 

S 

This  potential  satisfies  the  auxiliary  conditions  (3),  (4),  and  (5)  simply 
because  <?>s  does  and  because  all  these  conditions  are  linear.  Similarly  <p(P) 
satisfies  the  partial  differential  equation  (2)  at  all  points  of  the  half  space 
y  <  0  that  are  exterior  to  the  surface  S.  Thus  q>(P)  satisfies  equations  (2), 
(3),  (4),  and  (5)  regardless  of  the  nature  of  the  function  c(q).  This  func¬ 
tion  is  determined  in  such  a  way  that  cp(P)  satisfies  the  boundary  condition 
(6). 


Applying  the  boundary  condition  (6)  to  the  potential  of  (13)  requires 
evaluating  the  limits  of  the  spatial  derivatives  of  (13)  as  the  point  P 
approaches  a  point  p  on  the  surface  S  (figure  3)*  These  same  limits  are 
required  to  calculate  fluid  velocity  components  on  S.  Care  is  required 
because  the  derivatives  of  l/r,  where  r  =  r(P, q)  is  given  by  (9)>  become 
singular  as  the  surface  is  approached.  The  other  terms  that  comprise  <pg  in 
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equation  (8)  are  not  singular  as  P  -» p.  Thus  the  singularity  is  exactly  the 
same  as  that  encountered  in  references  [l]  and  [2],  and  the  results  of  these 
references  apply  here.  In  particular,  the  application  of  the  boundary  condi¬ 
tion  (6)  to  the  potential  (13)  gives 

-  2jt<Kp)  +  [<Ps(p><l)M<l)<iS  =  f(p)  (l4) 

S 

where  the  right  side  of  (6)  has  been  written  f(p)  to  show  its  dependence  on 
position  on  S.  Equation  (l4)  is  a  Fredholm  integral  equation  of  the  second 
kind  for  the  source  distribution  o(q).  Once  it  is  solved  for  <x(q)  the 
potential  at  any  point  is  given  by  (l3)  and  the  velocity  at  any  point  by 

v(P)  =  -  grad  <p(P)  =  -j^grad  [cpg(P,q)]a(q)  dS  (15) 

S 

The  evaluation  of  either  (l3)  or  (15)  at  a  point  P  =  p  on  S  requires  the 
limiting  process  discussed  in  references  [l]  and  [2). 

6.5  Logic  of  the  Numerical  Procedure 

It  should  be  noted  that  the  discussion  of  the  previous  section  applies  to 

the  problems  of  references  [l],  [2],  [3]  as  well  as  to  the  present  problem. 

In  fact  the  equations  of  section  6.4  apply  exactly  to  each  of  the  problems  of 

these  references  if  q>  is  set  equal  to  the  appropriate  point-source  function. 

s 

Similarly  the  method  of  numerically  solving  the  problem  defined  by  equation  (l4) 
and  subsequently  calculating  the  quantities  (13)  and  (15)  is  logically  identical 
for  all  problems.  This  method  is  described  in  great  detail  in  the  references. 

It  will  simply  be  outlined  here. 

The  body  surface  S  is  approximated  by  a  large  number  of  small  plane 
quadrilateral  surface  elements,  as  is  shown  schematically  in  figure  4.  The 
total  number  of  elements  used  to  approximate  the  surface  is  denoted  N,  and 
typical  elements  are  denoted  the  i-th  and  the  j-th,  as  shown  in  figure  4.  Sub¬ 
scripts  i  and  j  denote  quantities  associated  with  the  i-th  and  j-th  elements, 

•A 

respectively.  In  particular  n^  is  the  unit  normal  vector  to  the  i-th  element. 
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Figure  4.  -  Approximation  of  the  body  surface  by  plane  quadrilateral  elements. 


and  ASj  is  the  area  of  the  j-th. element.  On  each  element  the  source  density 
is  taken  as  constant.  The  value  of  source  density  on  the  j-th  element  is 
denoted  Oy  On  each  element  a  control  point  is  selected  where  the  integral 
equation  (l4)  is  required  to  be  satisfied  and  where  the  potential  and  velocity 
will  eventually  be  evaluated.  The  control  point  of  the  i-th  element  is  denoted 

li¬ 
lt  is  clear  from  the  above  that  the  potential  and  velocity  at  the  control 
point  of  the  i-th  element  due  to  the  j-th  element  are,  respectively, 

<t.  .a.  and  V.  .<7.  (l6) 

ij  j  ij  J 


where 


=  JJ ds 

(17) 

■  ■  iTerad  ds 

AS. 

0 
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t 


and  where  denotes  a  general  point  of  the  j-th  element  as  shown  in  figure  4. 
The  normal  velocity  at  the  control  point  of  the  i-th  element  due  to  all  the 
elements  is 


'ni 


N 

-  I  V 


(18) 


3=1 


where 


Aij  "  ni 


id 


(19) 


Application  of  the  normal  derivative  boundary  condition  at  the  control  points 
of  all  elements  then  gives 


N 

I  Vj  -  -  f(pt>’ 


i  «  1,  2,  ...  N 


(20) 


J=1 


where  f(pi)  is  the  right  side  of  (6)  or  (l4)  evaluated  at  p^.  (Recall  that 
d<p/dn  =  —  Vn. )  Equation  (20)  is  the  numerical  approximation  to  the  integral 
equation  (l4).  It  consists  of  a  set  of  linear  algebraic  equations  for  the 
values  of  source  density  on  the  elements.  Once  this  set  of  equations  is 
solved,  the  potentials  and  velocities  at  the  control  points  of  the  elements 
are  calculated  from 


N 

\  *  I  \iV 


i  =  1,  2,  ...  N 


d-i 

N 


(21) 


VI 


V3’ 


i  =  1,  2,  ...  n 


d-i 


All  of  the  above  is  quite  straightforward.  The  only  problem  is  the  cal¬ 
culation  of  the  integrals  in  (lY)  for  a  function  ©  given  by  (8).  Hat  is, 

s 

the  only  new  calculation  that  is  required  to  solve  the  problem  of  oscillatory 
motion  near  a  free  surface  is  the  integration  of  the  oscillating  point-source 
potential  and  its  derivatives  over  a  quadrilateral  element. 
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7.0  CALCULATION  OF  THE  OSCILLATING  POINT-SOURCE 
POTENTIAL  AND  ITS  DERIVATIVES 

7*1  General  Remarks 

As  stated  in  the  previous  section,  the  central  problem  of  this  method  of 
solution  is  the  integration  of  the  oscillating  point-source  potential  and  its 
derivatives  over  a  plane  quadrilateral  surface  element.  These  integrations 
must  be  done  numerically  and  they  must  be  done  very  rapidly,  because  a  total 
of  N^  such  integrations  are  required  to  produce  the  matrices  and  V„ , 

where  N  is  the  number  of  elements  used  to  approximate  the  body  surface  and 
is  thus  a  large  number.  The  problem  of  numerical  integration  divides  itself 
into  two  parts.  Hie  first  is  the  calculation  of  the  point-source  potential  and 
derivatives  for  general  values  of  the  parameters.  The  second  is  the  integra¬ 
tion  of  these  quantities  over  an  element  by  means  of  a  multipole  expansion. 

This  section,  7*0,  will  consider  the  first  of  these  parts,  and  section  8.0 
will  consider  the  second. 


7.2  Hie  Quantities  Calculated 

The  horizontal  distance  R  between  the  point  source  at  the  point  (a,b, c) 
and  the  field  point  (x,y, z)  is  given  by  equation  (ll).  The  vertical  distance 
between  the  field  point  and  the  image  point  (a, -b,c)  of  the  point  source  is 

h  =  -(y  +  b)  >  0  (22) 

(See  figure  2  for  an  illustration  of  these  quantities.)  Using  these  variables 
and  the  results  of  eqi  ation  (12)  the  oscillating  point-source  potential  (8) 
can  be  written 


*■ =  -  *{l ^e'^l)(vR)  (23) 

The  l/r  and  3/r^  terms  are  already  accounted  for  by  the  method  of  refer¬ 
ence  [l].  They  will  not  be  considered  here,  but  attention  is  directed  to  the 
last  two  terms.  It  should  be  noticed  that  x  and  z  do  not  enter  the 


» 
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above  equation  separately  but  occur  only  in  R,  i.e. ,  the  potential  is 
axi symmetric  about  a  vertical  line  through  the  point  source.  Thus,  only 
h-  and  R-derivatives  need  be  considered,  and  the  x-  and  z-derivatives 
are  easily  obtained  from  the  latter.  All  distances  in  the  last  two  terns 
of  (23)  enter  as  products  with  the  inverse  distance  v.  Such  products  are 
dimensionless  and  are  now  defined  as  new  variables.  In  particular  the 
dimensionless  horizontal  and  vertical  distances  are 

A  =  v  R 

(24) 

B  =  v  h 

These  are  essentially  distances  measured  in  wave  lengths  of  the  motion. 
Substituting  these  in  above,  together  with  t  =  v  u  gives  the  following 
form  for  the  oscillating  point-source  potential: 

%  -  7  +  ^  +  +  V  (25) 

where 


)-2/  ,-e  dt  - 

(26) 

0  Va2  +  (t  -  B)2 

%  =  <^(A,B)  =  2irie*Vl)(A) 

(27) 

This  form  of  cps  is  used  for  A  >  5-  For  A  <  5  the  terms  are  regrouped  in 
the  following  manner: 

<ps  =  7  +  ^  +  v[-<%  +  'PjJ  (28) 


where 


% 


CC 

=  VA’B>  ■ 2/ 


e-t  dt 


0  Va2  +  (t  -  B)2 


+  27re”BY  (A) 
o  ' 


(29) 


<Pj  =  q>j(A,B)  =  2irie‘BJo(A) 


(30) 
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The  value  A  =  5>  where  the  change  is  made  from  (25)  to  (28),  is  more  or 
less  arbitrary.  Form  (25)  is  more  convenient  at  large  values  of  A.  Form 
(28)  must  be  used  at  small  values  of  A,  because  and  <Pj  are  well- 

behaved  at  A  =  0,  while  and  q^  are  singular  there. 

It  is  the  potentials  qp^,  q^,  q^,  q>j  and  their  derivatives  that  are 
considered  in  this  report.  Only  the  A-derivatives  need  be  considered.  The 
B-derivatives  can  be  expressed  in  terms  of  the  potential  itself.  Specifically, 


=  -  q> 


q>  =  %  °r  <Pj 


_ . 

*  <P  + 


/a2  +  B2 


ffl  +  , 

Y  vr^  * 


9=  q^ 


The  additional  term  in  (32)  can  be  handled  by  the  method  of  reference  [l). 

For  the  integration  over  an  element  discussed  in  section  8.0  all  four 
of  the  above  potentials  must  be  considered.  In  the  present  section,  however, 
the  interest  is  in  formulas  for  expressing  the  point-source  potential  itself. 
The  two  potentials  qjg  and  q>j  each  consist  of  an  exponential  times  a 
Bessel  function,  and  these  functions  may  be  calculated  from  more  or  less 
standard  formulas.  (Reference  [6]  contains  a  rather  extensive  section  on 
Be s se 1  funct ions . ) 

Thus,  in  this  section  attention  is  concentrated  on  the  potentials  9^ 
and  q^  and  on  their  A-derivatives.  Generally  q^  is  of  interest  for 
A  >  5  and  9^  for  A  1  5,  but  both  have  been  investigated  over  a  wider 
range  of  A  so  that  the  cross-over  point  could  be  changed  if  desired.  The 
potentials  q^  and  q^.  and  their  A-derivatives  are  shown  in  figures  1,  8, 
9,  and  10,  and  numerical  values  are  given  in  tables  1  and  2.  The  nature  and 
general  magnitude  of  these  functions  are  illustrated  in  these  figures. 
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TABLE  2 

AND  ITS  DERIVATIVE  d<fy/dA 
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7. 3  Laguerre-Gauss  Quadrature  for  Large  Values 
of  Horizontal  Distance 


For  "large"  values  of  the  dimensionless  horizontal  distance  A  between 
the  point  source  and  the  field  point,  the  oscillating  point- source  potential 
is  used  in  the  form  (25),  and  the  potential  from  (26)  is  the  one  that 
must  be  evaluated.  Laguerre -Gauss  quadrature  is  used  to  evaluate  <p^  in  the 
form 


% 


OC 

-*f 


e-t  dt 


/  — ■  £ -  a  2  )  — . -d-  - 

0  Va  +  (t  -  B)2  j=l  Va2  +  (tj  -  B)J 


(33) 


The  constants  H.  are  the  weights  of  the  quadrature,  and  the  numbers  t . 

J  3 

are  the  abscissas  of  the  quadrature.  Both  Hj  and  tj  depend  only  on  M, 

the  order  of  the  quadrature.  For  the  present  application,  the  first  abscissa 

was  set  equal  to  zero,  i.e., 


tL  -  0 


W 


so  that 


2H 


<Pl 


•■J 


H 


1 


Va2  +  j=2  Va2  +  (t  -  b)2 

J 


2H 

vr 


f-J 


(35) 


H 


I 


j=2  Va2  +  (tj  -  B)2 


In  general  the  accuracy  of  Laguerre -Gauss  quadrature  improves  with  increasing 
M,  i.e.,  the  more  abscissas,  the  more  accuracy.  Ihe  term  with  abscissa  of 
zero,  i.e.,  the  first  term  of  (35)>  is  computationally  free,  because  it  may 
be  combined  with  the  l/r^-term  of  (25).  for  the  same  computational 

labor  the  form  (35)  should  be  more  accurate  than  the  form  of  Laguerre-Gauss 
quadrature  that  does  not  prescribe  any  abscissas,  because  the  former  uses  one 
more  abscissa.  This  plausible  hypothesis  was  verified  in  the  present  case  by 
a  large  amount  of  numerical  experimentation. 
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Reference  [7]  gives  the  abscissas  t^  and  weights  Hj  for  the  form  of 
Laguerre -Gauss  quadrature  that  has  zero  value  of  the  first  abscissa,  The 
results  for  M  =  1,  2,  3,  4,  5>  6  are  given  in  table  3«  These  were  used  in 
equation  (33)  or  (35)  and  the  results  were  compared  with  the  results  of  "brute 
force"  numerical  evaluations  of  «  and  its  A-derivative.  (The  approxima¬ 
tion  to  the  A-derivative  of  is  obtained  by  analytically  differentiating 
(33)  or  (35)*)  In  this  manner  computational  errors  were  obtained  for  various 
values  of  M  at  sets  of  values  of  A  and  B.  For  each  A  the  largest 
error  for  any  value  of  B  was  determined  and  designated  the  maximum  error 
at  that  A.  These  maximum  errors  are  presented  in  table  4  for  2  <  A  <  10. 
These  errors  should  be  compared  with  the  functions  <p^  and  which 

are  given  in  table  1.  Table  4  illustrates  the  very  rapid  decrease  of  error 
with  increasing  A  and  M.  If  an  allowable  error  is  defined  as  one  less 
than  0.001,  two  abscissas  are  sufficient  for  A  >6,  while  three  suffice  for 
A  >  4.  Thus,  in  the  present  context  A  becomes  "large"  at  a  value  of  four 
or  less.  The  smaller  the  value  of  A,  the  larger  is  the  value  of  M  that 
must  be  used  to  obtain  a  given  accuracy. 


For  the  present  purpose,  it  was  tentatively  decided  to  use  M  =  3  and 
to  restrict  use  of  this  formula  to  the  case  A  >  5.  For  larger  values  of  M 
the  computation  (35)  is  more  time-consuming  than  the  expansions  discussed 
later  in  this  section.  With  this  M  errors  in  and  dcp^/dA  are  bounded 
by  0.00023.  The  formula  for  is  explicitly 


% 


=  0.66666667  + _ 1.244oi69 _  _ 0.089316398 

Va2  +  B2  \Ja2  +  (1.2679492  -  B)2  \/a2  +  (4.7320508  -  B)2 


A  >  5 
(36) 


and  the  formula  for  dcj^/dA  is  obtained  by  differentiating  (36). 


7.4  The  Exponential-Integral  Expansion  for  Small  Values  of  the 
Ratio  of  Horizontal  to  Vertical  Distance 


For  small  values  of  the  dimensionless  horizontal  distance  A  between  the 
point  source  and  the  field  point,  the  oscillating  point-source  potential  is 
used  in  the  form  (28),  and  the  potential  cp^  from  (29)  is  the  one  that  must  be 
evaluated.  On  page  477  of  reference  [4]  two  forms  of  the  oscillating  point- 
source  potential  are  given.  By  comparing,  these  forms  it  is  evident  that  in 
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TABLE  3 


ABSCISSAS  AND  HEIGHTS  FOR  LAGUERRE-GAUSS  QUADRATURE  WITH 
ZERO  VALUE  OF  THE  FIRST  ABSCISSA 


5 

‘3 

H3 

M  =  I 

1.0 

0 

1.0 

M  =  2 

B9 

0 

ms^m 

mm 

2.0 

E^B 

M  =  3 

mm 

0 

0.33333333 

H 

1.26794919 

0.62200847 

m 1 

4.73205081 

0.04465820 

M  =  4 

0 

0.25 

1 

0.93582223 

0.62905268 

3.30540729 

0.11835639 

IMA 

7-75877048 

0.00259093 

M  =  5 

■a 

0 

0.20 

HI 

0.74329193 

0. 60120469 

B 

2.57163501 

0.18573233 

5.73117875 

0.01294285 

Bfl 

IO.9538943I 

_ 

0.00012013 

M  =  6 

1.0 

0 

0.16666667 

2.0 

0.61703085 

0.56401481 

3.0 

2.11296596 

0.23771357 

4.0 

4.61083315 

0.03056192 

5.0 

8. 39906697 

0.00103820 

6.0 

14.26010307 

0.00000484 
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MAXIMUM  ERRORS  IN  EVALUATING  THE  LINE  SOURCE  POTENTIAL  qfc  USING 
IAGUEERE - GAUS S  QUADRATURE  WITH  VARIOUS  NUMBERS  OF  ABSCISSAS 


the  present  notation  the  following  result  holds 


!HHK(k/v)B jo(v *)■*  (57) 

1  0 

where  the  designation  PV  signifies  that  the  principal  value  of  the  integral 
is  to  be  taken.  With  the  variable  change  k  =  vt  and  some  rearrangement 
(37)  yeilds  for  <j^ 


*  e'^J  (At) 


/e  u  \ 

-r* 


-a — r 

Jf?  +  B2  0 


e'BtJo(At)dt 


=  2FV 


n  e-BtJ  (At) 

/-rir- 


This  last  follows  from  the  fact  that  the  bracketed  term  in  (38)  is  zero  as  can 
be  verified  from  a  table  of  Laplace  transforms,  for  example,  page  1024  of 
reference  [6].  Hie  form  (39)  serves  as  the  basis  of  an  expansion  involving 
the  exponential  integral  Ei(B),  which  is  defined  as 


00  t 

Ei(B)  =  -  PV  J  ~ 


A  simple  variable  change  gives  the  result  that 


Q  =  PV  J  ^7-jr  dt  =  e-BEi(B) 


from  which  it  immediately  follows  that 


S  ■  (->"  ”f  ^  -  -  £r  ['"H 
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TO  expand  the  integral  in  (39) >  replace  Jq  by  its  power  series,  which  is 
convergent  for  all  values  of  the  argument.  This  gives 


% 


=  2 


y  /Af 

L  /k02\2 ; 

k=0  ''k*' 


f  t2ke-Bt 

J  1  - 1 


rv 


dt 


In  view  of  the  above  this  is 


From  (1*0) 


h  [Ei(B)3  -  r 


(45) 


w 


(45) 


Kepeated  use  of  (1*5)  in  (1*4)  gives  the  final  form  of  the  exponential-integral 
expansion.  Specifically, 


w  21c 

%  ■ 2 1  (i) 

k=0 


where 


CQ(B)  =  e_BEi(B) 


(46) 


(47) 


Ck(B) 


B 


2k 


22k(k.')2 


2k- 1 


-  £ 


7s=0 


bA+! 


(48) 


The  expansion  for  the  A-derivative  of  cp^  is  obtained  from  (4o)  by  differ¬ 
entiation.  The  bracketed  term  in  (48)  consists  of  the  difference  between  the 
quantity  e_BEi(B)  and  the  first  2k  terms  in  its  asymptotic  expansion  for 
large  B  (see  below). 


50 


2k 

The  factor  B  has  been  included  in  the  coefficient  functions  given  by 
(48)  to  minimize  their  variation  with  B.  That  this  is  the  correct  factor  to 
accomplish  this  may  be  inferred  from  the  limiting  forms  of  these  functions 
for  very  small  and  very  large  values  of  B.  For  large  B  the  asymptotic 
series  for  the  exponential  integral  gives 


ck0s) 


(2k  ~  1);  2k 

2  2V)2  B 


For  small  B  the  form 


B  — *  00 


m 


C  (B)  -»-■£§£ - ^  [1  +  0(B)]  ,  B  -*  0  (50) 

k  22k(k.')2 

is  obtained  directly.  Urns,  C^(B)  approaches  zero  as  l/B  for  large  B 

and  approaches  a  finite  value  for  small  B.  Division  of  these  quantities  by 
2k 

B  would  result  in  coefficient  functions  that  rapidly  approach  zero  for 
large  B  and  that  become  large  for  small  B.  If  the  limiting  forms  (49)  and 
(50)  are  used  in  the  series  (46),  in  both  cases  the  ratio  test  shows  converg¬ 
ence  for  A/B  <  1  and  divergence  for  A/B  >  1. 

Thus,  the  expansion  for  small  values  of  A  is  actually  an  expansion 
in  powers  of  A/B,  a: id  it  is  the  value  of  this  parameter  that  determines  how 
many  terms  of  the  expansion  must  be  used  to  obtain  a  certain  desired  accuracy. 
Clearly  small  values  of  A/B  may  correspond  to  rather  large  values  of  A  if 
B  is  sufficiently  large. 

The  approximate  form  of  (46)  that  is  used  for  computation  is,  of  course, 

%  =  2  I  <-l>VB>(t f  (5!) 

k=0 


Values  of  (fi  calculated  from  (51)  and  of  dqf^/dA  calculated  from  the 
derivative  of  (51)  have  been  compared  with  "brute  force"  numerical  evaluations 
for  various  orders  of  the  expansion,  i.e.,  various  values  of  M,  over  a  wide 
range  of  values  of  B  and  A/B.  For  each  value  of  A/B  the  maximum  error 
for  any  value  of  B  was  determined  and  designated  the  maximum  error  at  that 


value  of  A/B.  These  maximum  errors  are  shown  in  figure  11.  (The  functions 
themselves  are  shown  in  figures  9  and  10. )  The  errors  in  the  derivative 
have  been  multiplied  by  A  to  minimize  their  dependence  on  B  for  fixed 
A/B.  For  M  =  4  the  expansion  is  accurate  to  0.001  up  to  values  of  A/B 
of  rbout  0.4  to  0.5.  The  smallest  value  of  M  that  is  used  is  M  =  1. 

Notice  that  this  means  that  two  terms  of  the  expansion  (51)  are  used. 

If  subroutines  are  available  for  the  exponential  integral,  the  coef¬ 
ficient  functions  C^(B)  can  be  calculated  directly.  However,  such  sub¬ 
routines  are  rare.  To  fill  this  need,  the  functions  Ck(B)  have  been  fitted 
with  polynomials  in  the  ranges:  0  s  B  <  1,  1  <  B  <  2,  2  £  B  £  4,  4  <  B  <  8, 
and  8  <  B  <  10.  The  results  are  contained  in  Appendix  A.  All  the  Ck(B) 
have  individual  polynomial  fits,  except  C^B),  which  is  calculated  from 
the  relation 

C1(B)  =  £  [b2Cq(B)  -  1  -  B]  (52) 

which  xs  obvious  from  the  definition  (48).  This  method  of  calculating  C^B) 
not  only  saves  computing  time,  but  also  prevents  numerical  difficulty  in  the 
raultipole  expansion  at  small  values  of  R  or  A  (see  section  8.3).  The 
other  coefficient  functions  can  also  be  expressed  as  a  multiple  of  Cq(b) 
plus  a  polynomial,  but  the  expressions  contain  2k  +  2  terms  and  for  the 
larger  k  are  more  complicated  than  the  fitted  polynomials.  In  Appendix  A 

all  the  C  (b)  are  fitted  to  about  the  same  accuracy,  approximately  0.0001. 

2k 

When  these  functions  are  used,  they  are  multiplied  by  (A/b)  ,  where  A/B 

is  never  larger  than  one-half.  Thus,  less  accuracy  is  required  for  the 
functions  with  larger  values  of  k,  and  further  investigation  will  very 
probably  show  that  one  or  two  terms  may  be  omitted  from  the  polynomial  fits 
for  ck(B),  k  =  2,  3>  4.  This  change  will  reduce  computing  time. 

For  B  >  10,  the  coefficient  functions  are  calculated  using  the 
asymptotic  expansion  of  e-BEi(B)  for  large  B.  This  expansion  is 


B  — *  00 


(53) 


» 
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The  accuracy  of  an  asymptotic  expansion  is  usually  optimized  if  the  expansion 
is  terminated  with  the  smallest  term.  For  B  =  10,  the  tenth  and  eleventh 
terms  of  (53)  are  the  smallest,  but  the  accuracy  is  insensitive  to  term 
number  between  six  and  ten.  For  the  present  purposes  six  terms  are  used  for 
0^,  Cg,  and  C y  and  eight  terms  are  used  for  C^.  As  mentioned  above, 
use  of  the  asymptotic  expansion  leads  to  seme  cancellation  in  the  bracketed 
term  in  .(48).  The  resulting  expressions  for  the  coefficient  functions  are:. 


co<B> 


120 

T 


C1(B)  »  £  [B\(B)  -  1  -  b] 
°2(B)  ■&[•*£] 


c3(b)  *  0 


(54) 


c4(b)  *  o 

The  approximations  (54)  are  less  accurate  for  the  higher  values  of  k,  but 
the  lessening  of  accuracy  with  increasing  k  is  consistent  with  the  decreasing 
accuracy  requirement  that  is  discussed  above. 


7«5  Expansions  for  Finite  Values  of  the  Eatio  of  Vertical 
to  Horizontal  Distance 

The  expansion  of  the  previous  section  calculates  <f^(A,B)  accurately  for 
all  values  of  the  parameter  B/A  larger  than  some  value  in  the  range  2.0  to 
2.5,  say  B/A  =  2.5  for  definiteness.  What  now  remains  to  be  devised  is  a 
method  of  calculating  cp^(A,  B)  for  the  range  0  <  B/A  <  2.5*  It  turns  out 
that  several  expansions  are  required  to  do  this,  each  of  which  is  based  on  a 
particular  value  of  B/A. 
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It  is  desired  to  expand  <fy(A,B)  about  a  particular  value  of  B/A, 
vtoich  is  denoted  p,  i.e.,  the  expansion  is  about  B  =  pA.  All  the  expan 
sions  of  this  type  are  based  on  the  formula 


%(A,pA)  +  2 


(B/A)-p 


Au 

e  du 


^1  +  (u  + 


0  Vl  +  fu  +  p  y 


(55) 


which  is  obtained  from  (29)  by  straight-forward  manipulation.  First  def 


me 


=  B  -  pA 


(56) 


The  function  1/  Jl  +  (u  +  p)2  i 
the  form 


is  expanded  as  a  power  series  about  u  =  0  in 


==  ■  X  ak(p),,‘ 


\Jl  +  (u  +  p)2  k=0 

where  the  first  few  a^p)  are  as  follows 


(57) 


\A  +  P2  aQ(p)  =  1 


>A  +  P2  a  (p)  - - 

1  +  P 


l 


1  + 


P2  a2(p)  = 


P2  -  1/2 
(I  +  P2)2 


'll  +  p2  a,(p)  =  -  £ — Z-0/2)p 
5  (1  +  P2)5 

JL  +  P2  at  (p)  =  p ■  —  ?P  +,?/8 

(1  +  p2)4 


v/l  +  p2  a_(p)  =  -  ^  +  (!5/8)p 

5  (1  +  P2)5 


(58) 
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Now  (55)  becomes 


%(A,B) 


e 


-€ 


%(A,pA)  +  2 


k  Au, 
u  e  du 


(59) 


The  integral  in  (59)  is  a  standard  form  with  value 


*4 


k  Au,  e _ 

‘  e  dU  =  TkTT 
0  A  a=o 


\  ,  k:  k-A  .  (-l)1^:  /  e  .v 

L  (_1)  xr^-xy.»e  +Ajrfr”(e  -1) 


(60) 


This  gives  for 


q^(A,B)  =  e 


-€ 


£!2£!£,,,> 


_k-A 


\A  e 

i_i  (k  —  A)-1 

k=l  A  A=0 


00 

2(e*  - 1)  £  <-Dk  -fir  \(p) 


k=0 


(61) 


This  can  be  rearranged  to  give  the  final  result 

M 

<H,(A,B)  .  e'^tA.pA)  +  Y  ('1>\(p>flt+2(e,(i)  (&) 

k=0 


where 
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A  +  P2  bL(p)  =  — ^ 

1  +  p 


A  +  P2  b  (p)  =  i  ^ 

2  5  (1  +  P2)2 


(65) 
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L  .  2  .  ,  _  v  l  8o^  -  24o2  + 
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L  .  2  .  /.x  1  8p5  -  40o5  + 

'1  +  p  b5(p)  =  5T  (1  + 


and  where  the  functions  f^(e)  are  defined  as 


fk(c)  -  I  <-)J  Ik^TT  «J 


k-2 

-  t-1’"'1  *‘£  - 1  (-p*  £ 

6  L  *=°  . 


i  —  —  + 

k  kuT* 


ihe  derivative  of  is  given  by 


(A,B)  «  e~€  )  pi 


W-l 

»%(*,(*)  + |f  k(A,e»)]Uj  ^  (-i)\(P)fW(£)  (f)  (fi5) 

'  1.  A 


where 


'1  +  P  c0(p)  =  -  2p 


(l  +  p-  c  (p)  =  ■  -  g 

1+P 


fl  +  P  c2(p)  =  3 


(TTTr 


(66) 


(1  ♦  P2)3 


^ +  p2  c^p)  =  ? 

(i  +  p2) 


h  p 

8p  -  12  p  +  1 


Notice  that  the  infinite  series  that  would  be  required  to  give  equality  in  (62) 
and  (65)  have  been  replaced  by  finite  series  for  practical  computation.  The 
upper  limit  M  of  the  summation  for  the  potential  is  the  order  of  the  expan¬ 
sion.  A  total  of  M  +  1  terms  involving  powers  of  e  are  added  to  the  basic 
term  in  the  potential  expansion  (60),  and  M  +  2  terms  in  powers  of  €  sure 
added  in  the  derivative  expansion  (63).  (Both  expansions  are  terminated  with 
the  same  fk(e).) 

If  the  expansions  (62)  and  (65)  are  considered  infinite  series,  the 
question  of  their  convergence  arises.  From  analysis  of  special  cases  and 
general  order  of  magnitude  arguments,  it  appears  that  the  series  convergence 
for  | e/A  |  <  1. 


The  above  formulas  permit  expansions  up  to  M  =  U  to  be  written  down 
for  any  expansion  parameter  p.  (in  fact  M  =  5  is  included  for  the  potential 
expansion. )  One  value  of  p  is  special,  namely  p  =  0.  For  this  value, 
every  other  term  in  the  expansion  is  zero,  and  for  the  same  computational 
effort  roughly  twice  as  large  a  value  of  M  can  be  used  for  p  =  0  as  for 
other  values  of  p.  Expansions  have  been  worked  out  for  three  values  of 
p:  p  =  0,  p  =  1,  p  =  7/b.  These  expansions  are  designated  the  near-zero, 
near-one,  and  near-seven-quarters  expansions,  respectively.  They  appear  to 
cover  the  range  0  s  B/A  <  2.3  adequately,  but  possibly  a  few  small  ranges 
require  new  expansions.  For  these  values  of  p,  the  first  few  coefficients 
in  the  expansions  are  given  in  table  5«  Specifically,  table  5  contains  bk(l), 
and  ck(7A)  for  k  =  0,  1,  2,  3,  ^  5  and  bk(0)  and 
0^(0)  for  k  =  1,  2,  3>  . ..,  9*  With  these  coefficients  the  near-zero  expan¬ 
sion  can  be  written  down  with  up  to  five  nonzero  terms  in  the  summations  —  the 
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same  number  as  the  other  expansions.  For  k  >  5  the  general  formulas  for 
b^(p)  and  c^(p)  are  not  given  above,  because  these  terms  are  used  only  for 
p  =  0.  For  nonzero  p  it  is  more  efficient  computationally  to  derive  more 
expansions  than  to  extend  existing  expansions  to  hi$i  values  of  k. 

Hie  accuracy  of  the  near-zero  and  near-one  expansions  for  various  orders 
M  has  been  determined  by  comparing  values  of  and  8<p^/8A  calculated  by 
the  expansions  with  "brute  force"  numerical  evaluations  for  a  considerable 
range  of  values  of  A  and  B/A.  For  each  value  of  B/A  the  largest  error 
that  occurs  for  any  A  has  been  determined  and  designated  the  maximum  error 
at  that  value  of  B/A.  These  maximum  errors  are  shown  in  figures  12  and  13 
for  the  near-zero  and  near-one  expansions,  respectively.  The  errors  in  the 
derivative  have  been  multiplied  by  A  to  minimize  their  dependence  on  A  for 
fixed  B/A.  An  exhaustive  error  study  has  not  been  made  for  the  near-seven- 
quarters  expansion.  However,  from  the  form  of  the  general  expansion,  it  is 
felt  that  the  errors  for  this  expansion  are  nearly  the  same  as  those  for  the 
near -one  expansion  if  the  two  expansions  are  compared  at  equal  values  of 
e/A  =  B/A  —  p.  That  is,  figure  13  can  serve  approximately  as  the  error  plot 
for  the  near-seven-quarters  expansion,  if  the  abscissa  scale  is  translated  to 
make  7/4  lie  where  1  does  now. 

For  small  k  the  functions  f^(e)  are  evaluated  from  the  middle  line 
of  (64),  i.e., 

fL(e)  =  e”e 

ef2(e)  =  1  -  e"e  (67) 

e2f5(e)  =  2[e"e  +  €  -  l] 

It  will  be  noticed  that  the  functions  f^(e)  enter  the  expansions  (62)  and 
(65)  multiplied  by  the  powers  of  e  given  in  (67).  For  higher  values  of  k, 
this  method  of  evaluation  is  cumbersome.  Accordingly,  the  functions  f^.(e ) 
for  k  =  3,  4,  5  and  6  have  been  fitted  with  polynomials  over  a  range 
-2.5  <  €  <  2.5.  The  functions  fk(e)  for  k  =  8,  and  10,  have  been 
fitted  over  the  range  0  <  e  <  2. 5*  These  latter  functions  are  used  only 
in  the  near-zero  expansion  for  which  e  cannot  be  negative.  The  polynomial 
fits  of  these  functions  are  given  in  Appendix  B. 


To  evaluate  expansions  (62)  and  (65)  explicit  formulas  for  the  leading 
terms  are  required.  Specifically,  the  functions 

G(A,p)  =  cp^pA) 
and 

K(A,p)  =  pq^(A,  pA)  +  ~  Pq^(A,pA)J  (68) 

are  needed  in  the  range  0  <  A  <5*  These  two  functions  have  been  fitted  with 
polynomials  in  the  range  1  <  A  <  5  for  p  =  0,  1,  and  7/4.  The  results 

are  given  in  Appendix  C.  For  the  range  0  <  A  <  1  improved  polynomial  fits 
were  obtained  by  writing  the  above  functions  in  the  forms 

G(A,p)  =  g(A,p)  +  2e-pAJo(A)ln  ( |  A j 

H(A,p)  =  h(A,p)  -  2e"pAJ1(A)ln  (|  a)  (69) 

+  2e  ^JjA)  J 

The  functions  g(A, p)  and  h(A, p)  are  not  singular  for  small  A,  and  it  is 
these  functions  that  are  fitted  with  polynomials  in  the  range  0  <  A  <1.  The 
results  are  also  given  in  Appendix  C. 

It  can  be  seen  from  equation. (69)  that  these  expansions  are  singular  for 
A‘=  0.  With  p  fixed  the  condition  A  -»  0  implies  that  also  B  -» 0  and 
thus  that  the  source,  its  image,  and  the  field  point  approach  each  other  and 
approach  the  free  surface.  Clearly,  singularity  is  expected  for  this  situation. 
Both  these  expansions  and  the  multipole  expansion  of  the  following  section  fail 
for  this  condition. 

If  the  expansions  are  used  for  values  of  A  as  large  as  five,  the 
ranges  of  the  fits  in  Appendix  B  are  appropriate  as  long  as  |  e/A |  is  no 
greater  than  one-half.  If  it  should  prove  necessary  to  use  values  of  |g/aI 
greater  than  one-half,  either  the  range  of  the  polynomial  fits  must  be  expanded, 
which  is  quite  straightforward,  or  the  expansions  must  be  restricted  to  values 
of  A  less  than  some  value  less  than  five.  In  the  latter  case,  the  Laguerre- 
Gauss  formulas  must  be  used  down  to  the  new  limit  on  A,  and  this  can  be 
arranged  simply  as  can  be  seen  in  table  4.  However,  it  is  desirable  to  keep 
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values  of  U/a|  as  low  as  possible  so  that  relatively  few  terms  of  the 
expansions  (62)  and  (65)  are  needed  for  good  accuracy. 

7.6  A  Sample  Range  Criterion  for  the  Expansions 

The  formulas  of  the  previous  parts  of  this  section  can  be  used  to  calculate 
and  its  A-derivative  to  any  degree  of  accuracy.  The  greater  the  accuracy 
that  is  required,  the  greater  is  the  number  of  expansions  that  are  needed, 
i.e.,  the  greater  is  the  number  of  values  of  p  that  must  be  used  with  (55)* 

The  accuracy  requirement  basically  determines  idiich  expansion  and  how  many 
terms  of  it  are  to  be  used  at  a  particular  set  of  values  of  A  and  B.  How¬ 
ever,  the  translation  of  this  requirement  into  an  actual  logical  decision 
routine  for  selecting  the  expansion  and  term  number  is  a  nontrivial  task.  To 
illustrate  this  procedure  an  example  is  presented  here. 

For  illustrative  •  oses  it  is  assumed  that  it  is  required  to  calculate 

cp^  and  d<f^/dA  with  an  error  not  exceeding  0.001  in  absolute  value.  (Compare 
the  values  of  these  functions  given  in  table  2  and  figures  9  and  10. )  This 
seems  to  be  a  reasonable  error  criterion  for  applications.  Possibly,  it  is  too 
conservative.  (A  decision  routine  for  a  different  error  criterion  could  be 
worked  out  from  the  error  curves  of  figures  11,  12,  and  13 . )  It  is  further 
assumed  that  only  the  expansions  discussed  in  the  previous  parts  of  this 
section  are  to  be  used  with  only  the  terms  that  have  been  presented,  The 

calculation  is  to  cover  the  complete  range  of  B  and  all  A  up  to  five.  In 

all  cases  the  fewest  possible  terms  are  to  be  used  to  conserve  computing  time. 

The  error  curves  shown  in  figures  11,  12,  and  13  are  used,  but  the 
results  are  adjusted,  always  in  the  conservative  direction,  to  simplify  the 
decision  routine.  The  results  are  shown  in  figure  l4,  whose  abscissa  is  A 
and  whose  ordinate  is  B/A.  The  coordinate  plane  is  divided  into  regions, 
in  each  of  which  a  certain  number  of  terms  of  a  particular  expansion  is  used 
to  calculate  and  Scp^/dA.  In  all  cases  the  number  of  terms  is  specified 

by  M,  where  this  quantity  has  the  same  meaning  as  it  does  in  the  definitions 
of  the  expansions  that  are  given  in  the  previous  parts  of  this  section.  The 
error  curves  of  the  expansions  have  been  conservatively  simplified  to  make 
the  boundaries  of  the  regions  in  figure  l4  straight  lines,  whose  equations 
are  given.  It  can  be  seen  that  the  only  regions  that  are  not  handled 
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satisfactorily  are  three  small  regions  having  small  values  of  A  near  B/A 
values  of  0.6,  1.4,  and  2.1,  respectively.  Hxese  regions  could  be  handled 
by  means  of  one -term  expansions  about  those  values  of  p.  Possibly  an 
increase  of  the  error  criterion  would  cause  these  regions  to  disappear. 


8.0  INTEGRATION  OVER  AN  ELEMENT  USING  A  MULTIPOLE  EXPANSION 


8.1  Geometric  Quantities  Associated  with  an  Element 

To  carry  out  the  method  of  solution  outlined  in  section  6. 5,  it  is 
necessary  to  develop  formulas  for  calculating  the  potential  and  velocity  at 
points  in  space  due  to  a  unit  source  density  on  one  of  the  plane  quadrilateral 
elements  used  to  approximate  the  body  surface  (see  figure  4).  This  requires 
integrating  the  formulas  for  the  oscillating  point  source  potential  and  its 
gradient  over  a  general  quadrilateral  element.  Certain  geometric  quantities 
associated  with  the  quadrilateral  are  used  to  perform  this  -integration. 

These  quantities  are  defined  in  this  section. 

Figure  5  shows  a  general  quadrilateral  element.  The  coordinates  of  the 
centroid  of  the  area  of  the  quadrilateral  are  xq,  yQ,  Zq.  Integrations  are 
performed  in  a  coordinate  system  based  on  the  element.  5he  variables  of  this 
coordinate  system  are  5,  q,  where  the  and  q-axes  lie  in  the  plane  of 
the  element  and  the  £-axis  is  normal  to  this  plane.  The  origin  of  this  coordi¬ 
nate  system  is  taken  as  the  centroid,  and  the  unit  vectors  along  the  axe3  of 
the  system  are: 


^-(x0t-y0,z0) 

h 

0 

FREE  SURFACE 

1  (*.  y.  i) 

X 

<^AS 

R° 

*o) 

(a.b.c) 

Figure  5.  -  Illustration  of  quantities  used  to  obtain  the  multipole  expansion. 


(70) 


5-axis: 

alli  +  a12d  +  ai3k 

^  A  ^ 

q-axis : 

a21z  +  a22J  +  a23k 

5 -axis: 

hi1  +  V  +  a33k 

where  i,  j,  k  are  the  unit  vectors  along  the  x,  y,  z  axes,  respectively. 
A  point  with  coordinates  x,  y,  z  has  coordinates  q,  £  in  the  system 
based  on  an  element  where  these  coordinates  are  related  by 


6  -  au<*  -  *„)  +  V<y  "  yo>  +  a13(z  ‘  zo) 


-1  .  a21(x  -  xo)  ♦  a22(y  -  yo)  +  a2J(z  -  zj 


(71) 


£  =  a3i^x  - xJ +  a^(y "  yJ  +  a^z  -  zJ 
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The  maximum  dimension  of  the  quadrilateral  is  denoted  t.  Its  value 
determines  the  size  of  the  quadrilateral.  Qhe  shape  of  the  quadrilateral  is 
expressed  by  the  normalized  moments  of  its  area.  Specifically,  the  normalized 
moments  are: 


where  the  integral  is  over  the  area  of  the  element.  The  order  of  a 

particular  moment  is  defined  as  the  sum  of  its  subscripts  rc+n.  In  general, 

the  size  of  the  moments  decreases  with  increasing  order.  The  zeroth  order 

moment  I0Q  is  the  normalized  area  of  the  element,  and  the  second  order 

moments,  I0_,  I,.,  and  I00,  are  the  normalized  "moments  of  inertia".  The 
20  Ij.  02 

first  order  moments,  1^  a.nd  1^^,  are  zero  because  the  origin  of  coordinates 
is  the  centroid  of  the  area. 
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8.2  The  Multipole  Expansion 


It  will  be  recalled  that  the  potential  at  a  point  (x,y,z)  due  to  a 
unit  oscillating  point  source  at  a  point  (a,b,c)  is  of  the  form 

<Ps  =  r  +  ^+V<P  (73) 

In  integrating  this  expression  over  an  element  the  two  l/r  terms  are  treated 
by  the  method  of  reference  1.  The  term  of  interest  here  is  of  one  of  two 
forms.  Either 


vcp’(A,B)  =  v[  B)  +  <^(A,B)]  (74) 

or 

vq>'(A,B)  =  v[-q^(A,B)  +  9j(A,B,]  (75) 


where  the  individual  potentials  qp^,  q^,  cp^,  and  q)j  are  given  by  equations 
(2 6),  (27),  (29),  and  (30),  respectively.  In  this  section 


9  =  <p(a,b) 


(76) 


is  used  to  denote  any  one  of  the  four  individual  potentials  of  (74)  and  (75) 
which  must  be  integrated  over  an  element.  The  quantities  A  and  B  depend  on 
the  coordinates  a, b,c  of  the  point  source.  Specifically, 


A  = 


v/(x  -  a)2  +  (z  -  c)2  = 


vR 


B  =  -v(y  +  b)  =  vh 


(77) 


The  process  used  to  integrate  the  quantity  V9  over  an  element  is  known  as 
the  multipole  expansion.  In  general  terms  it  proceeds  as  follows.  The  point 
(a, b,c)  where  the  source  is  located  is  taken  as  a  general  point  on  the  quad¬ 
rilateral  element  as  shown  in  figure  5*  This  point  is  expressed  in  the  element 
coordinate  system  by  using  a,b,c  in  place  of  x,y, z  in  equation  (7l)» 

Because  the  point  is  on  the  element  the  resulting  f  is  zero.  By  this  means 
the  quantities  A  and  B,  and  thus  the  quantity  V9,  are  expressed  in  terms 
of  the  coordinates  x  ,  y  ,  z  of  the  centroid  of  the  element  and  the 
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coordinates  £.  q  of  a  point  on  the  element  in  its  own  coordinate  system. 

This  expression  for  v<p  is  expanded  in  a  two-variable  Taylor  series  in  powers 
of  |  and  q.  The  resulting  series  consists  of  certain  derivatives  of  <p 
evaluated  at  the  centroid,  i.e.,  at  the  point  |  =  q  =  0,  multiplied  by  powers 
of  |  and  q.  The  derivatives  of  9  are  constants  with  regard  to  integration 
over  the  element,  and  integration  of  the  powers  of  £  and  q  gives  the  nor¬ 
malized  moments,  equation  (72). 


The  subscript  0  is  used  to  denote  quantities  evaluated  at  the  centroid, 
e.g (ifcp/dA )q,  etc.  In  particular 


Aa  =  v  \/(x  —  x  )^  +  (z  —  z  =  vR 
0  0  0  o 


B0  =  “V(y  +  yo}  =  Vho 


(78) 


Here  Rq  and  hQ  are,  respectively,  the  horizontal  and  Arertical  distances 
between  the  field  point  (x,  y,  z)  and  the  image  point  (xq, -yq,Zq)  of  the 
centroid  of  the  element. 


In  carrying  out  the  details  of  the  raultipole  expansion  it  turns  out  to  be 
more  efficient  to  expand  the  velocity  components  rather  than  the  potential. 

Let  the  integrated  potential  of  the  element  be 


4>  = 


V  JJ 

£3 


(79) 


where  9  may  be  any  of  the  potentials  cp^,  9^,  9^, 
components  associated  with  this  potential  are 


or  9j- 


The  velocity 


V  = 
x 


V  = 
y 


v 


z 


(80) 


b6 


These  components  are  integrated  over  the  element  in  the 
manner.  The  potential  can  be  expressed  in  terms  of  the 
Specifically, 


above-described 
y-  component  V  . 

y 


-**  jit- 

a  r*  -L 


if 


q>=  ^ 


or 


<t>  =  - 


-  V 

v  y 


(8l) 


if 


9  =  %  or  <p 


The  integral  term  of  (8l)  can  be  absorbed  with  the  integral  of  the  l/r^  term 
of  (73)  and  evaluated  by  the  method  of  reference  1. 

In  deriving  explicit  formulas  for  the  velocity  components  (80)  use  is 


made  of  the  fact  that  the  first  order  area  moments, 


Z10  and 


01’ 


are  zero. 


Use  is  also  made  of  the  fact  that  9  is  an  axisymmetric  solution  of  Laplace's 
equation  in  cylindrical  coordinates,  i.e.,  that 


(82) 


The  manipulations  required  to  derive  the  formulas  of  the  multipole  expansion 
are  rather  lengthy  and  are  not  included  here.  Instead  only  the  results  are 
presented. 
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■  •  v*a*»Cj4| 


The  expressions  for  V  and  V  are  identical  to  the  above  with  superscripts 

y  “ 

x  replaced  by  superscripts  y  and  z,  respectively.  Before  explicit  formulas 
for  the  various  •'uantities  in  (83)  are  presented,  some  discussion  of  the  gen¬ 
eral  form  of  the  expansion  seems  in  order.  The  expansion  (83)  is  the  sum  of 
two  expansions  in  powers  of  (vt):  the  "b  array"  whose  rows  are  multiplied  by 
successive  B-derivatives  of  <p,  and  the  "c  array"  whose  rows  are  multiplied 
by  successive  B-derivatives  of  (8<p/8A).  All  rows  of  both  arrays  are  infinite 
series  in  powers  of  (t/RQ),  except  the  first  row  of  the  "b  array",  which  con¬ 
sists  of  a  single  term.  The  complete  expansion  is  thus  in  powers  of  the  two 
variables  (vt)  and  (t/R  ).  The  coefficients  b  and  c  (with  any 
superscript  x,y,  or  z)  have  subscripts  that  denote  the  powers  of  those 
variables  that  the  coefficient  multiplies.  Specifically,  the  first  subscript  m 


48 


denotes  the  power  of  (vt),  and  the  second  subscript  n  denotes  the  power 

of  (t/R  ).  The  coefficients  b  and  c  are  each  linear  combinations  of 
'  o  mn  mn 

the  normalized  moments  of  order  equal  to  two  less  than  the  sum  of  the  subscripts 

of  the  coefficient.  (That  is,  the  sum  of  the  subscripts  of  the  coefficient  is 

two  greater  than  the  sum  of  the  subscripts  of  the  normalized  moments  from  which 

it  is  formed. )  Thus,  each  coefficient  b  or  c  is  designated  as  being 

mn  mn 

of  a  particular  order,  which  is  simply  the  order  of  the  normalized  moments 
used  to  form  it.  The  chief  decrease  in  size  of  the  coefficients  in  (83)  is 
with  increasing  order.  Note  that  coefficients  of  a  fixed  order  are  all  on 
the  same  "antidiagonal"  of  the  arrays  in  (83).  This  is  similar  to  the  multipole 
expansion  described  in  reference  [3]. 

In  particular,  the  terms  of  (83)  involving  b2Q  and  c^q  are  the  zeroth 
order  terms  of  the  expansion.  These  coefficients  depend  on  the  normalized 
area  Iqq  of  the  element.  The  first  order  coefficients  b^Q,  c^,  and  c^ 
depend  on  the  first  order  moments  and  I  .  Thus,  these  coefficients 

are  zero  if,  as  in  the  present  case,  the  expansion  is  about  the  centroid  of 
the  element.  They  have  been  included  in  (83)  to  make  clear  the  form  of  the 
expansion.  The  coefficients,  b^,  b^,  c^,  c^,  and  c22  are  the  second 
order  coefficients.  The  expansion  is  used  in  two  forms.  Either  the  zeroth 
order  terms  are  the  only  ones  retained,  or  the  zeroth  order  terms  and  all 
second  order  terms  are  retained.  In  the  latter  case,  expansion  (83)  is  used 
as  written.  Higher-order  terms  are  not  included  because  the  complexity  of  the 
terms  increases  rapidly  with  order. 

First  define  the  auxiliary  quantities 


x  —  x 


a  = 


R 


P  = 


y  -  y„ 


z  —  z 


R 


7  = 


R 


s-l” 

o 


a12^  * 


5  =  R-  "  V 

O 


(84) 


(85) 
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G1  =  *12*20  +  a22Ill  ' 


al2  ll  +  a22I02 


H  «  al2Gl 


a22°2 


j-*j„  +  e-  i^-  i. 


20 


02 


Then  the  coefficients  in  expansion  (35)  are  as  given  in  the  tables  below: 


m 

i 

b(x) 

mn 

b(y) 

mn 

b(z) 

mn 

! 

0 

0 

I00 

0 

0 

0 

0 

0 

H 
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D 

0 

"l 
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*1 

m 

n 

CW 

mn 

c(y) 

mn 

c(z> 

mn 

2 

0 

aIoo 

0 

rIoo 

2 

1 

0 

0 

0 

2 

2 

?b(x) 

2b31 

0 

2b(l) 

°3l 

3 

0 

0 

0 

0 

3 

1 

(2b4o^  “  anGi  “  a21G2) 

-  \  (J  “  2J2) 

Ho’  -  a13Gl 

a25G2) 

■ 

0 

-J1 

/40 

The  many  zeros  in  these  tables  greatly  simplify  the  computation. 


The  B-derivatives  in  (83)  can  be  simply  related  to  the  quantity  that  is 
differentiated.  It  is  evident  from  definitions  (27)  and  (30)  or  directly 
from  relation  (31)  that 


if 


<P  =  %  or  <Pj  (91) 


For  <p,  and  the  B-derivatives  can  be  computed  successively  from  the 
definitions  (2 6)  and  (29)  or  more  directly  from  relation  (32).  !lhe  results 
are  as  follows: 


if  cp  = 


or 


% 


&  =  + 

dB2 


<£2  = 


dB' 


<p- 


<p  - 


<p  - 
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2B 


Vrl  (vr,)5 


vr 


2B  _ 2_ 


6bc 


o  y 

5bSa 


1  (vr,)5  (vr,)' 
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(vr,)' 


M  A  2A 
“  (vr,)5 


d5cp 

dB2dA 


dg>  +  2A  +  6ab 
^  (vr,)5  (vr,)5 


(92) 


where  as  before 


Vr,  =  /a2  + 


B 


(93) 


Expressions  similar  to  (92)  could  be  derived  for  the  higher  derivatives,  sc 
that  in  principle  expansion  (83)  can  be  carried  to  any  desired  order.  The 
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derivatives  in  (92)  are  all  that  are  needed  if  expansion  (83)  is  terminated 
with  second  order  terms. 

The  formulas  of  this  section  permit  the  potential  and  velocity  at  a 
general  point  in  space  due  to  a  general  quadrilateral  element  to  be  evaluated 
in  terms  of  the  potential  and  A-derivative  of  the  potential  due  to  a  point- 
source  at  the  centroid  of  the  quadrilateral.  These  two  quantities  are 
evaluated  by  the  formulas  of  section  7-0.  Thus,  the  multipole  expansion 
reduces  the  quadrilateral  source  to  a  point  source  as  far  as  computation  is 
concerned  and  is  effective  for  all  range  of  A  and  B. 

8.3  Behavior  of  the  Multipole  Expansion  at  Small  Values 
of  Horizontal  Distance 

From  the  form  of  expansion  (83)  it  appears  that  the  expansion  becomes 
singular  as  RQ/t  approaches  zero,  i.e.,  as  the  horizontal  distance  between 
the  centroid  of  the  element  and  the  field  point  becomes  small.  However,  it 
is  clear  that  if  the  potential  <p  is  an  analytic  function,  as  are  all  9 
considered  for  the  present  purpose,  this  singularity  must  be  apparent  and  not 
real.  It  is  possible  that  the  apparent  singularity  might  lead  to  numerical 
difficulty,  and  calculational  procedures  must  be  designed  to  avoid  such  prob¬ 
lems.  It  turns  out  that  to  ■'avoid  any  singularity  the  multipole  expansion 
must  include  either  all  terms  of  a  particular  order  or  no  such  terms.  The 
procedure  adopted  in  this  report  is  such  that  the  apparent  singularity  causes 
no  difficulty.  Of  course,  expansion  (83)  cannot  be  used  if  RQ/t  is  exactly 
zero.  Thus,  a  test  is  made  and  special  limiting  formulas  are  used  if 
|RQ/t|  <  e,  where  e  is  taken  as,  say,  10~\  The  discussion  below  shows 
that  the  expansion  formula  may  be  used  for  small  values  of  RQ/t  just  greater 
than  any  nonzero  e. 

The  potentially  troublesome  terms  of  (83)  are  those  containing  negative 
powers  of  RQ/t.  For  the  present  discussion  it  suffices  to  consider  only 
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those  terms  and  to  ignore  the  others.  The  sum  of  the  singular  terms  of  (83) 
are  (dropping  the  subscripts  o) 


Vising.)  =-(^f)(vt)3 

4X)  (if  w 

+  (at?)  (vt>3  (l) 

For  finite  values  of  v  and  B  =  vh  taking  the  limit  as  R  -»  0  means  that 
A  =  vR  and  thus  A/B  are  both  small.  Thus,  the  potential  is  considered  only 
in  the  form  (28),  i.e.,  <f^  and  <Pj  are  used  in  expansion  (83),  not  and 


the  exponential-integral  expansion  using  an  older  of  the  expansion  M  =  1 
(figure  1*0.  Thus,  equation  (51 )  gives  for  t'le  potential  and  its  derivative 


^  *  2  [C0(E)  -  cl(e)  (£)2].  o(rI‘) 

sr  -  [-  5  (|)]+  °(r3> 


(95) 


where  the  O(R^)  and  0(R3)  terms  in  (95)  merely  illustrate  the  magnitude 
of  neglected  terms,  since  only  the  terms  in  square  brackets  are  retained  for 
M  =  1.  The  coefficient  C^(b)  is  not  calculated  independently  but  is  obtained 
from  Cq(b)  by  equation  (52).  Thus,  the  above  expressions  are 


%  ■  - 1  tvR>2] +  Hr  (if + 

J  (96) 

^!.-Co(B)(vR)  +  4-i(£)*0(R3) 
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,;n<  *«»•*■ 


and 


If  (96)  sind  (92)  are  used  in  (9*0,  the  result  is  (dropping  the  0(R  ) 
0(lc  )  terms  that  contribute  nothing  in  the  limit) 


Vising. )  =  - 
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Vx(sing.)  = 
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(98) 


(x) 

31 


The  two  terms  in  curly  brackets  in  (98)  will  be  computed  separately  in  practice. 
Both  terms  have  singularities  of  order  i/R  with  coefficients  that  are  the 
terms  in  square  brackets.  These  singularities,  which  are  shown  below  to 
cancel  analytically,  must  also  cancel  numerically.  Since  the  singularity  is 
only  first  order,  this  can  be  arranged  with  moderate  caution  in  the  calcula- 
tional  procedure.  This  has  been  verified  explicitly  by  numerical  experimenta¬ 
tion  at  R  /t  =  e  =  10-^. 
o' 

There  are  two  kinds  of  terms  in  (98).  The  terms  that  do  not  involve 
Co(B)  are  exact  except  for  round-off  caused  by  the  finite  word  length  of 
the  computer.  The  terms  containing  Cq(b)  have  an  additional,  much  larger 
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error  due  to  the  polynomial  fit  of  Cq(b).  To  show  these  two  types  of  terms 
explicitly,  (98)  is  written  in  the  form 


* 

"  - 

Vising.)  =  (vt)5 

Co(B) 
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-kd5i  f  c22 
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(vt)  c^' 
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B  +  1 
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2 

(vr^ 


+  0(R) 


(99) 


From  the  table  of  coefficient  functions  above,  it  can  be  seen  that 


(100) 


and  the  same  statement  is  true  if  the  superscript  x  is  replaced  by  y  or  z. 
Thus,  the  coefficient  of  (t/R)CQ(B)  in  (99)  vanishes,  and  thus  errors  due  to 
fitting  Cq(b)  are  not  critical.  This  is  entirely  due  to  the  fact  that  (B) 
is  computed  in  terms  of  Cq(b)  from  (52).  If  C^(B)  were  computed  separ¬ 
ately,  Cq(b)  and  C^(b)  would  have  independent  errors  and  the  required 
cancellation  would  not  occur.  There  would  then  be  numerical  difficulty  for 
small  R.  To  show  the  disappearance  in  the  limit  of  the  second  term  in  (99) 
it  should  be  noticed  that 


1 

vr. 


y/k2  +  B2 


=  |  +  0(R2) 


(101) 


This,  together  with  (100)  gives  the  desired  result.  Thus,  in  the  limit  as 
R->  0 


Vising.)-*  (vt)  c 


(x)  fR\  _  B  +  5 

31  [VB>  b2  J 


(102) 


with  similar  expressions  for  V  and  These  expressions  may  be  used  for 

!  R  /t|  <  e. 

I  o' 
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Methods  of  computing  the  potential  (f>j,  which  is  given  by  (30),  are  not 
discussed  in  section  7*0,  because  it  is  felt  that  standard  expressions  will 
be  adequate.  In  particular,  for  moderate -to-small  values  of  A,  it  seems 
reasonable  to  use  the  rapidly  convergent  Taylor  series  for  Jq(A)  in  powers 
of  A.  In  this  case  the  apparent  singularities  in  the  multipole  expansion 
can  be  shown  to  cancel  analytically  and  numerical  difficulty  at  small  R 
should  be  no  worse  than  for  cp^. 

8.4  Accuracy  of  the  Multipole  Expansion 

The  accuracy  of  the  multipole  expansion  was  investigated  by  using  it  to 
calculate  a  large  number  of  cases  and  comparing  the  results  with  the  results 
of  "brut*  force"  numerical  integrations.  In  formulating  error  curves  one 
problem  was  to  be  conservative  without  being  too  conservative.  The  procedure 
used  to  generate  these  error  curves  is  described  below. 

The  potential  was  considered  in  the  form  f28),  i.e.,  divided  into 
and  <Pj.  (The  l/r  terms  were  not  considered. )  Thus  the  potential  was 
divided  into  its  real  part  -v%  and  its  imaginary  part  v<Pj/i.  In  the 
derivative  expressions  (9l)  and  (92)  that  enter  into  the  expansion  (83)  exact 
values  of  cp  and  (dcp/dA)  were  used,  so  that  any  errors  are  those  a r  i 

with  the  multipole  expansion.  The  accuracy  of  the  multipole  expansion 
on  the  parameters  (hQ/t),  (RQ/t)  and  (vt)  and  also  on  the  shape  of  the 
element  and  on  the  direction  of  the  field  point  with  respect  to  the  centroid 
of  the  element.  In  calculating  error  curves,  the  maximum  error  with  respect  to 
element  shape,  direction,  and  velocity  component  was  determined  as  a  function 
of  the  three  parameters  above.  Specifically,  a  set  of  values  of  (hQ/t), 

(RQ/t),  and  (vt)  was  selected.  Calculations  were  performed  for  six  differ¬ 
ent  elements  and  for  a  variety  of  directions.  The  largest  error  in  any  velocity 
component  for  a  ny  of  the  elements  and  directions  was  determined  by  inspection. 
This  was  done  for  the  potentials  vcf^  and  vq)j  separately  to  obtain  a 
"maximum  real  part"  and  a  "maximum  imaginary  part"  error.  The  square  root 
of  the  sum  of  the  squares  of  these  two  errors  is  the  quantity  eventually 
plotted.  It  has  the  nature  of  a  "maximum  absolute  value"  error,  and  it  is 
conservative  because  the  "maximum  real  part"  and  "maximum  imaginary  part" 
errors  seldom  occur  for  the  same  velocity  component,  the  same  element,  and 
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the  same  direction.  The  above  procedure  was  repeated  for  ranges  of  values  of 
(hQ/t),  (R Q/t),  and  (vt)  to  obtain  the  maximum  error  curves  presented  here. 

Maximum  error  curves  were  calculated  for  the  case  t&ere  all  second-order 
multipole  terms  were  retained,  i.e.,  where  expansion  (83)  is  useu  as  written, 
and  also  for  the  case  \diere  only  zeroth-order  multipole  terms  are  retained, 
i.e.,  only  the  b^  and  c terms  of  (83).  Figure  15  shows  error  curves 
for  the  second-order  multipole  expansion,  and  figure  l6  shows  error  curves  for 
the  zeroth-order  multipole  expansion.  Values  of  (vt)  up  to  three  were  con¬ 
sidered.  For  (vt)  =  3  the  maximum  element  dimension  is  approximately  equal 
to  one-half  the  wave  length  of  the  motion.  No  higher  values  of  (vt)  were 
considered,  because  it  was  felt  that  at  such  frequencies  the  assumption  of  a 
constant  source  density  on  the  element  would  break  down.  If,  however,  this 
assumption  remains  valid,  higher  frequencies  can  be  considered  by  subdividing 
elements  to  reduce  the  maximum  value  of  (vt)  below  three  (see  next  section). 
A  comparison  of  figures  15  and  l6  for  (hQ/ t )  =  2  and  (hQ/ t )  =  3  shows  the 
very  large  improvement  In  accuracy  that  is  obtained  by  using  the  second- 
order  multipole  expansion  instead  of  the  zeroth  order.  This  fact  was  also  a 
factor  in  the  decision  not  to  proceed  to  higher  orders  for  the  multipole 
expansion. 

If  a  definite  error  criterion  is  adopted,  figures  15  and  l6  permit  ranges 
of  validity  to  be  established  for  the  zeroth  and  second-order  multipole  expan¬ 
sions.  For  definiteness,  it  is  assumed  here  that  an  absolute  error  criterion 
of  0.001  has  been  adopted.  With  this  criterion  the  ranges  may  be  defined: 

a.  For  (hQ/t )  >  4,  zeroth-order  multipole  is  sufficient  for  all 
(RQ/t)  and  all  (vt)  <  3«  (As  can  be  seen  from  figure  l6c, 
possibly  (hQ/t)  >  4.2  is  required  for  small  (RQ/t),  but  this 
is  a  detail. ) 

b.  For  4  >  (hQ/t)  >  2,  second-order  multipole  gives  the  required 
accuracy  for  all  (RQ/t)  and  all  (vt)  £  3* 

c.  For  (hQ/t)  <  2  second-order  raultipole  is  not  sufficiently  accurate 
for  all  (RQ/t)  and  all  (vt)  £  3»  However,  the  accuracy  improves 
very  rapidly  with  decreasing  (vt)  and  with  increasing  (RQ/t). 
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In  the  range  (hQ/t)  <  2  the  integration  over  an  element  is  accomplished  by- 
subdividing  the  element  into  subelements,  each  of  -which  has  a  smaller  value 
of  t.  Thus,  when  calculating  the  effect  of  the  subelement,  the  effective 
value  of  (vt)  is  reduced  compared  to  that  of  the  original  element.  In  most 
cases  the  values  of  (h^/t.)  and  (R ^/t)  are  increased  relative  to  those  of 
the  original  element.  As  can  be  seen  from  figure  15,  all  three  of  these 
changes  in  parameter  values  lead  to  a  reduction  of  the  error  for  the  second- 
order  multipole.  The  computation  of  the  effect  on  itself  of  an  element  adja¬ 
cent  to  the  free  surface  may  require  special  handling.  This  is  acceptable 
from  the  standpoint  of  computing  time,  because  such  cases  represent  a  small 
fraction  of  the  total  computation.  However,  the  use  of  a  second-order  multipole 
with  element  subdivision  will  give  accurate  results  for  the  important  case  of 
the  effect  of  one  element  adjacent  to  the  free  surface  on  another  such  element. 
It  is  only  for  elements  near  the  free  surface  that  subdivision  is  required, 
and  thus  it  is  used  rather  infrequently  in  the  computation  scheme. 

8.5  A  Scheme  for  Subdividing  an  Element 

The  subdivision  of  an  element  for  use  with  the  multipole  expansion  at 
values  of  (hQ/t)  less  than  two  may  be  accomplished  in  a  variety  of  ways. 

From  a  computational  standpoint  the  process  of  subdividing  an  element  con¬ 
sists  of  the  calculation  of  the  sets  of  geometrical  quantities  that  define  the 
subelements.  Specifically,  the  following  quantities  must  be  calculated  for 

each  subelement:  the  normalized  moments  I  of  the  area,  the  coordinates 

mn 

of  the  centroid,  the  maximum  dimension  t,  and  the  components  a^,  a^,  etc. 
of  the  unit  vectors  along  the  axes  of  the  subelements '  coordinate  systems. 

(in  reference  1  the  set  of  these  a  's  is  called  the  transformation  matrix.) 

If  further  subdivision  is  required  to  reduce  values  of  (vt)  to  acceptable 
values,  i.e. ,  if  the  subelements  must  themselves  be  subdivided,  the  coordinates 

systems  must  also  be  calculated.  An  efficient  subdivision  scheme  is  one 
that:  (l)  obtains  ail  the  above  quantities  with  as  little  computation  as 
possible,  (2)  reduces  values  of  (vt)  as  much  as  possible,  and  (5)  can  be 
iterated  in  a  straightforward  manner  to  subdivide  the  subelements.  Many 
schemes  are  possible.  One  scheme  that  possesses  these  requirements  is  out¬ 
lined  here. 


\ ,  k  =  1,2, 5,^  of  the  corners  of  the  subelements  in  their  own  coordinate 
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A  typical  element  to  be  subdivided  is  shewn  in  figure  6.  The 
coordinates  of  its  four  corner  points  in  the  element  coordinate  system  are 
k  =  1»2, 3 A-  It  is  assumed  that  the  maximum  dimension  t  of  the 
element  is  taken  as  the  longer  of  its  two  diagonals  and  that  the  |-axis 
of  the  element  coordinate  system  is  taken  parallel  to  this  longer  diagonal. 
Thus,  if  (as  shewn  in  figure  6)  the  longer  diagonal  is  between  the  points 
(lrnL)  and  then  =  tj and  t  ®  ({^  -  The  origin  of 

the  element  coordinate  system  is  taken  as  the  centroid  of  the  element. 

(All  of  the  above  is  consistent  with  the  geometrical  considerations  of 
reference  1.  The  only  required  change  is  making  the  |-axis  parallel  to 
the  longer  diagonal.) 


Conceptually,  the  subdivision  scheme  consists  of  bisecting  each  side  of 
the  quadrilateral  with  a  point  and  drawing  a  line  connecting  each  of  these 
four  points  with  the  midpoint  of  the  longer  diagonal.  As  shown  in  figure  6, 
this  process  yields  four  subelements,  which  are  labelled  1,  2,  3,  4  to 
denote  the  corner  point  they  contain.  Subelements  1  and  3  have  the  same  shape 


Figure  6.  -  Subdivision  of  an  element  by  the  use  of  midpoints. 
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and  orientation  as  the  original  element,  and  each  is  reduced  by  one -half  in 
linear  dimension.  Subelements  2  and  k  are  parallelograms.  Each  subelement 
has  one  diagonal  parallel  to  the  longer  diagonal  of  the  original  element. 

The  5-axes  of  the  coordinate  systems  of  all  subelements  are  thus  taken 
parallel  to  the  5-axis  of  the  coordinate  system  of  the  original  element. 

Thus,  the  components  a^,  a^,  etc.,  of  the  unit  vectors  along  the  axes 
of  the  coordinate  systems  of  all  subelements  are  identical  with  the  same 
quantities  for  the  original  element,  and  no  additional  computation  is  required 
for  these  quantities.  Subelem_nts  1  and  3  have  the  same  normalized  moments 
as  the  original  element  and  have  values  of  t  exactly  half  that  of  the  original, 
element.  Normalized  moments  must  be  computed  for  the  parallelogram  subelements 
2  and  4,  but  in  subsequent  subdivision  of  the  parallelograms  the  normalized 
moments  are  identical  for  all  later  subelements.  The  length  of  the  longer 
diagonal  t  of  one  of  the  parallelogram  subelements  cannot  be  predicted  in 
advance.  The  diagonal  parallel  to  the  5-axis  is  exactly  half  the  value  of  t 
for  the  original  element,  but  the  other  diagonal  may  be  longer.  In  the  worst 
case  the  other  diagonal  of  a  parallelogram  subelements  may  be  only  slightly 
shorter  than  the  value  of  t  for  the  original  element.  However,  in  subsequent 
subdivisions,  the  values  of  t  for  the  subelements  of  the  parallelogram  are 
exactly  half  the  value  of  t  for  the  parallelogram. 

Explicit  formulas  for  the  geometric  quantities  associated  with  the  sub¬ 
elements  1,  2,  3,  and  4  of  figure  6  show  how  little  computation  is  required 
by  the  subdivision  process.  First  define 


“p  “  I  ITT"  * 
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e2  "  5.-6, 
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Then  the  geometric  quantities  associated  with  the  subelements  in  figure  6  are 
as  follows: 
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Finally,  from  the  method  of  subdivision  it  is  clear  that  the  coordinates  of 
the  corner  points  of  the  subelements  of  figure  6  are  simple  averages  of  the 
coordinates  of  the  corner  points  of  the  original  element. 
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EXPRESSIONS  FOR  EVALUATING  THE  COEFFICIENT  FUNCTIONS  IN  THE 
EXPONENTIAL-INTEGRAL  EXPANSION  FOR  VALUES  OF  THE  ARGUMENT 

FROM  ZERO  TO  TEN 
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c0(B) 


e"B  InB 

+  0.57723459 
+  0.42188715  B 

-  0.45420044  B2 
+  0.18830916  B3 

-  0.036074946  B4 
+  E(B) 


|E(B)j  <  25  *  10' 


E(l)  =  -  0.000019 


UBS2 


cq(b)  =  0.69717880 


+  0.30182470  (B  -  1) 

-  0.63786604  (B  -  l)2 
+  0.48345948  (B  -  l)5 

-  0.22434376  (B  -  l)k 
+  0.050233747  (B  -  l)5 
+  E(B) 


|E(B)|  <  20  *  10' 

E(l)  =  0.000004 | 
E(2)  =  0.000004 


2  S  B  i  4 


Co(B)  =  0.67048886 

-  0.17123280  (B  -  2) 

-  0.034727768  (B  -  2)2 
+  0.042573269  (B  -  2)5 

-  0.014448463  (B  -  2)h 
+  0.0019073418  (B  -  2)5 
+  E(B) 


I  e(b)|  <  25  •  io" 

E(2)  =  0.000006  | 
E(l)  =  0.000006 


4  s  B  s  8 


co(B) 


0.35953823 

-  0.10948767  (B  -  4) 

+  0.023664245  (b  -  4)2 

-  0.0031370348  (B  -  4)5 
+  O.OOOI8870981  (B  -  4)4 
+  E(B) 


|  e(b)|  <  30  •  io"6 

E(4)  =  -  0.000014 
E(8)  =  -  o.ooooi4 


...  ix  1*4***. a 


8  §  B  5  10 


Cq(b)  =  0. 14772609 

-  0.022641305  (B  -  8) 

+  0.0033149587  (B  -  8) 2 

-  0.00032011296  (B  -  8)5 
+  e(b) 


|e(b)|  <  15  •  io"6 

e(8)  =  -  0.000005 
E(10)  =  -  0.000005 


Cjte)  =  [B2e_BEi(B)  -  1  -  B] 


In  all  ranges  of  B  C^B)  is  calculated  from  the  identity 
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2  £  B  £  4  C5(B)  =  -  0.092^91866 

-  0.028463143  (B  -  2) 

+  0.00033224846  (B  -  2)2 
+  0.0074613867  (B  -  2)3 

-  0.0011826565  (b  -  2)^ 

+  e(b) 

4  §  B  £  8  C?(B)  =  -  0.10729389 

+  0.023850510  (B  -  4) 

+  0.018257392  (B  -  4)2 

-  0.0039381087  (B  -  4)5 

-  0.000014366499  (B  -  4)^ 
+  0.000036616026  (B  -  4)5 
+  E(b) 

8  £  B  £  10  C5(b)  =  0.061968114 

+  0.023687746  (b  -  8) 

-  0.0086533615  (B  -  8)2 
+  0.00080471344  (B  -  8)5 
+  E(B) 


|  E(B)|  <  20  *  lO-6 


E(2)  =  -  0.000005 
E(4)  =  0.000004 


|E(B)|  <  40  *  lO-6 

IE(4)  =  0.000030 
E(8)  =  0.000028 


I  e(b)|  <  25  •  10"6 

E(8)  =  -  0.000009 
E(10)  =  -  0.000004 


C4(B)  =  [B®e‘^Ei(B)  -  B7  -  B6  -  2B5  -  6B^  -  24b5  -  120B2  -  702B  -  5040] 


0  £  B  £  1 


1  £  B  £  2 


CU(B)  =  -  0.034191950  1 

1  e(b)|  <  20  •  io‘6 

-  0.0047257587  B 

-  0.0011700246  B2 

J  E(l)  =  0.000014 

+  E(B) 

Cj^B)  =  -  0.040121243  1  E(b)|  <  25  •  lo"6 

-  0.0069648968  (b  -  l) 

'  -  0.0026223112  (B  -  l)2 

E(l)  =  -  0.000019 

+  E(B) 

E(2)  =  +  0.000017 

67 


2  §  B  S  4 


4  s  B  5  8 


8  §  B  £  10 


c4(b)  =  -  0.049763691 

-  0.0117^7816  (B  -  2) 

-  0.005035^531  (B  -  2)2 
+  0.001603^521  (B  -  2)3 
+  E(B) 


I  E(B)|  <  50  *  10*6 


E(2)  =  -  0.000038 
E(4)  =  -  0.000032 


Ck(B)  =  -  0.080550890 

-  0.011972214  (B  -  4) 

+  0.0048073235  (b  -  4)2 
+  0.0027910827  (B  -  4)3 

-  0.00072846533  (B  -  4)^ 

+  0.000046982416  (b  -  4)5 
+  e(b) 

C^(b)  =  -  0.011289474 

+  0.034544199  (B  -  8) 

-  0.0023406480  (B  -  8)2 

-  0.00057919267  (B  -  8)5 

+  E(B) 


I  e(b)|  <  20  •  io“6 

E(4)  =  -  0.000009 
E(8)  =  -  0.000005 


|  e(b)|  <  30  •  io-6 


E(8)  =  -  0.000024 

E(io)  =  -  0.000021 


68 


APPENDIX  B 


POLYNOMIAL  FITS  FOR  THE  COEFFICIENT  FUNCTIONS  APPEARING  IN  THE 
EXPANSIONS  FOR  FINITE  VALUES  OF  THE  RATIO  OF  VERTICAL 
TO  HORIZONTAL  DISTANCE 
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APPENDIX  C 

POLYNOMIAL  FITS  FOR  THE  LEADING  TERMS  IN  THE  EXPANSIONS  FOB  FINITE 
VALUES  OF  THE  RATIO  OF  VERTICAL  TO  HORIZONTAL  DISTANCE 
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Near-Zero  Expansion  (Derivative) 
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1  a  S  2.1  H(A,0) 


2.1  a  «  3  h(a,0) 
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Near-One  Expansion  (Potential) 
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Figure  9.  -  The  potential  V 


to  o 


Figure  14.  -  Ranges  of  values  of  A  and  B/A  where  certain  numbers  of  terms  of  the  several 
expansions  are  to  be  used  to  obtain  an  accuracy  of  0.001  in  computing  the 
potential  and  its  A-derivative. 


ERROR  IN  VELOCITY 


.0024 


90 
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Figure  15.  -  Continued,  (c)  ho/t  =  2. 
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